rm(list=ls())

# set working directory
setwd("")

library(foreign)
library(readstata13)
library(ggplot2)
library(devtools) 
install_github("susanathey/causalTree") 
library(causalTree)
library(gbm)

data <- read.dta13("dta_files/ind.dta")

#Causal trees (Athey)

# Bias scale
tree <- causalTree(biasscale ~ female + agedeciles + education + partnership + parent  + hhmembers,
                   data = data, treatment = data$treat,
                   split.Rule = "CT", cv.option = "CT", split.Honest = T, 
                   cv.Honest = T, split.Bucket = T, 
                   xval = 5, cp = 0, minsize = 20, propensity = 0.5)
opcp <- tree$cptable[,1][which.min(tree$cptable[,])]
opfit <- prune(tree, opcp)
opfit
rpart.plot(opfit)
dev.copy2pdf(file="output/Fig7a.pdf", width = 7, height = 5, title = "Behavioral bias")

# Right-wing scale
tree <- causalTree(rwscale ~ female + agedeciles + edu + partnership + parent  + hhmembers,
                   data = data, treatment = data$treat,
                   split.Rule = "CT", cv.option = "CT", split.Honest = T, 
                   cv.Honest = T, split.Bucket = T, 
                   xval = 5, cp = 0, minsize = 20, propensity = 0.5)
opcp <- tree$cptable[,1][which.min(tree$cptable[,])]
opfit <- prune(tree, opcp)
opfit
rpart.plot(opfit)
dev.copy2pdf(file="output/Fig7b.pdf", width = 7, height = 5)

# Populism scale
tree <- causalTree(popscale ~ female + agedeciles + educat + partnership + parent  + hhmembers,
                   data = data, treatment = data$treat,
                   split.Rule = "CT", cv.option = "CT", split.Honest = T, 
                   cv.Honest = T, split.Bucket = T, 
                   xval = 5, cp = 0, minsize = 20, propensity = 0.5)
opcp <- tree$cptable[,1][which.min(tree$cptable[,])]
opfit <- prune(tree, opcp)
opfit
rpart.plot(opfit)
dev.copy2pdf(file="output/Fig7c.pdf", width = 7, height = 5)

# Refugee rejection scale
tree <- causalTree(refscale ~ female + agedeciles + educat + partnership + parent  + hhmembers,
                   data = data, treatment = data$treat,
                   split.Rule = "CT", cv.option = "CT", split.Honest = T, 
                   cv.Honest = T, split.Bucket = T, 
                   xval = 5, cp = 0, minsize = 20, propensity = 0.5)
opcp <- tree$cptable[,1][which.min(tree$cptable[,])]
opfit <- prune(tree, opcp)
opfit
rpart.plot(opfit)
dev.copy2pdf(file="output/Fig7d.pdf", width = 7, height = 5)



